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ABSTRACT  Tensioned  or  compliant  structures,  such  as  fabric  shelter  systems 
and  ocean  cable  systems,  are  composed  of  very  flexible  structural  components  that 
exhibit  a  high  degree  of  geometrical  nonlinearity  and  analytical  complexity.  Structural 
analysis  of  these  systems  is  a  two- phase  procedure:  (1)  the  determination  of  the  static 
equilibrium  configuration  or  prestressed  state  of  the  structure;  and  (2)  the  determination 
of  structural  response  relative  to  the  equilibrium  configuration  due  to  static  and  dy¬ 
namic  in-service  loads.  This  report  addresses  static  solution  algorithms  for  reliable 
phase-one  solution  procedures.  In  contrast  to  the  standarad  Newton-Rapbson  (NR) 
incremental/iterative  algorithm,  the  automated  dynamic  relaxation  (ADR)  solution  al¬ 
gorithm  is  a  pure  iterative  algorithm  and  it  is  also  simple  to  implement.  The  ADR 
algorithm's  performance  is  globally  convergent  in  test  problems  having  multiple  solu¬ 
tions,  whereas  NR  algorithms  are  known  to  be  only  locally  convergent.  The  ADR 
algorithm  can  be  successfully  utilized  to  form  the  basis  of  structural  analysis  software 
for  specialized,  highly  nonlinear  problems  involving  tensioned  structures. 
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INTRODUCTION 


The  problem  addressed  in  this  research  is  the  lack  of  a  robust  com¬ 
putational  procedure  for  determining  the  initial  equilibrium  configura¬ 
tion  and  prestress  state  of  tensioned  structures.  The  solution  of  this 
problem  is  widely  regarded  as  a  deterrent  to  efficient  engineering  design 
and  analysis  of  tensioned  structures  in  the  Navy  (i.e.,  ocean  cable  and 
membrane  structures  and  land-based  expeditionary  and  architectural/ten¬ 
sioned  fabric  structures).  Several  authors  have  discussed  and  documented 
this  technical  problem;  they  include  Webster  (1977,  1984),  Liu  (1977), 
Peyrot  and  Goulois  (1979),  Shields  and  Zueck  (1984),  and  Shugar  (1987). 
Compliant  structures  such  as  these  generally  may  not  be  handled  by  linear 
analysis  methods,  and  successful  nonlinear  analysis  methods  usually  require 
modified  solution  algorithms . 

A  promising  algorithm  for  application  to  structures  that  exhibit 
strong  nonlinear  structural  behavior  was  addressed  in  this  investigation. 

It  is  known  as  the  automated  dynamic  relaxation  (ADR)  algorithm.  The 
ADR  algorithm  was  compared  to  existing  solution  algorithms  for  tensioned 
structures.  It  possesses  some  attractive  features  which  provide  for 
monitoring  and  control  of  the  stability  of  the  solution  process.  Often 
the  solution  process  exhibits  pathological  behavior  when  a  tensioned 
structure  approaches  a  slack  configuration. 

Existing  solution  algorithms  based  upon  the  Newton -Raphson  method 
have  a  subjective  nature.  They  require  that  a  guess  be  made  for  the 
pretension  in  the  structure  to  begin  the  solution  procedure.  The 
robustness  of  these  methods  is  adversely  affected  by  this  requirement 
because  convergence  to  the  correct  solution  sometimes  depends  on  the 
accuracy  of  the  guess.  This  runs  contrary  to  the  goal  of  a  foolproof 
algorithm. 
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A  set  of  cable  test  problems  was  designed  to  numerically  evaluate 
the  robustness  of  the  algorithms  studied.  In  these  problems,  the  start¬ 
ing  conditions  were  designed  purposely  to  be  onerous  to  test  the  ability 
of  the  algorithms  to  seek  the  correct  static  equilibrium  position  when 
starting  from  rather  arbitrary  configurations.  In  some  cases,  multiple 
static  equilibrium  solutions  exist,  and  the  goal  is  to  find  the  global 
solution. 

Objective 

The  objective  of  this  report  is  to  present  the  automated  dynamic 
relaxation  solution  algorithms  for  nonlinear  static  problems  in  a  manner 
suitable  for  implementation  into  structural  analysis  software.  This 
report  is  largely  based  on  a  previous  paper  (Shugar,  1988)  but  also 
includes  additional  background  information,  discussion  of  different  ADR 
algorithms,  and  additional  numerical  comparison  studies  with  general 
purpose  ocean  cable  analysis  computer  programs. 

Background 

The  technology  of  computational  methods  in  ocean  structural 
engineering  was  reviewed  by  Shugar  and  Armand  (1986)  in  which  solution 
methods  for  geometrically  nonlinear  ocean  cable  structures  were  dis¬ 
cussed.  In  particular,  the  solution  of  static  mooring  system  problems 
has  largely  motivated  this  study  of  the  dynamic  relaxation  solution 
method.  Hence,  a  lengthy  discussion  of  several  existing  numerical 
solution  procedures  and  mooring  computer  programs  is  included  in  the 
appendix.  The  computer  programs  discussed  include: 

MOORING  (NAVFAC  DM-26.5,  1985) 

FLEETMOR  (Palo  and  Karnoski,  1986) 

STATMOOR  (Cox,  1982) 

SOUPLE  (Peyrot,  1980) 

SEADYN  (Webster,  1976) 

SEASTAR  (Pawsey  and  Nour-Omid,  1988) 
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Three  classes  of  computer  programs  are  defined  as  follows:  (1) 
programs  limited  to  standard  Navy  ship  mooring  practice,  (2)  programs 
suited  to  special  purpose  mooring  problems,  and  (3)  general  purpose 
finite  element  programs  aimed  at  a  wide  class  of  structural  analysis 
problems . 

The  last  two  types  of  programs  listed  above  are  more  closely 
related  than  the  first  in  that  they  are  based  on  the  direct  stiffness 
method  from  the  theory  of  matrix  analysis  of  structures,  which  provides 
them  with  the  potential  for  general  purpose  capability.  The  first  type 
is  very  specialized  and  limited  to  ship  mooring  problems  in  shallow 
water  where  mooring  line  hydrodynamic  drag  force  is  insignificant.  With 
respect  to  the  static  mooring  problem,  however,  all  three  types  of  pro¬ 
grams  have  the  common  requirement  of  solving  a  nonlinear  system  of 
static  equilibrium  equations. 

The  Newton  iteration  solution  algorithm  is  most  commonly  employed 
to  solve  the  nonlinear  system  of  equations  in  the  computer  programs 
discussed.  None  of  these  programs  employ  the  dynamic  relaxation  solu¬ 
tion  algorithm,  although  it  has  been  used  successfully,  mostly  in  the 
U.K. ,  for  general  purpose  static  analysis  of  very  flexible  or  compliant 
structures  including  tensioned  fabric  structures  and  offshore  riser  and 
cable  structures.  Mathematical  and  engineering  descriptions  of  the 
dynamic  relaxation  method  are  given  by  Crisfield  (1986),  although  the 
context  is  restricted  to  the  solution  of  linear  equations. 


DERIVATION  OF  DYNAMIC  RELAXATION  INTEGRATION  FORMULAS 

It  is  presumed  that  a  continuous  cable  structural  system  has  been 
spatially  discretized  using  standard  simple  straight  truss  or  curved 
isoparametric  finite  elements.  Then,  the  system  of  N  static  nonlinear 
algebraic  equations  to  be  solved  for  the  unknown,  N-dimensional  displace¬ 
ment  vector  x  is: 

F,(x)  =  F  (1) 

~i~  ~e 
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where  F,(x)  is  the  internal  force  vector  which  is  a  function  of 

structural  system  configuration  x,  and  F@  is  the  applied,  static 

external  force  vector  acting  on  the  structural  system. 

The  dynamic  relaxation  method  begins  by  converting  the  static 

problem  (Equation  1)  into  a  structural  dynamics  problem  using  a  psuedo 

time  variable,  t,  and  a  step  load  form  for  the  static  load  vector, 

F  (t).  Thus, 

~e 


M  x  (t)  +  C  x  (t)  +  FA  [x  (t)]  =  Fe  (t) 


(2) 


Here  M  and  C  are,  respectively,  artificial  mass  and  artificial  damping 
matrices  which  are  defined  arbitrarily  as: 


M  = 

p  I 

r  M 

(3a) 

c  = 

C  I 

(3b) 

where  I  is  the  NxN  identity  matrix,  and  the  scalars  p  and  c  are  artifi¬ 
cial  mass  and  artificial  damping  parameters,  respectively.  These  param¬ 
eters  will  subsequently  be  chosen  to  control  convergence  of  the  dynamic 
relaxation  iterative  procedure.  The  artificial  inertia  and  artificial 
damping  forces  are  stabilizing  terms  for  the  left-hand  side  of 
Equation  2.  When  calculation  of  the  internal  force  becomes  numerically 
ill-conditioned,  these  terms  serve  to  temporarily  regularize  the 
calculation.  The  psuedo  equation  of  motion  (Equation  2)  is  to  be 
integrated  for  x(t).  After  sufficient  time  has  elapsed,  the  dynamic 
solution  x(t)  is  expected  to  converge  to  the  static  solution  sought 
which  is  denoted  by  x. 

Temporal  integration  is  accomplished  by  discretizing  Equation  2 
using  the  following  central  difference  approximation  formulas, 


x(tn-l/2) 


•n-!/2  . 


(4a) 
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and, 


x(V  5n  =  K  (jn+1/2  “  2n’1/2) 


(4b) 


where  h  is  the  time  step  size.  Also,  averaging  the  discrete  velocities 
in  Equation  4a  time  steps  n+1/2  and  n-1/2  yields: 


i(t  )  xn  =  |  (in+1/2  -  in_1/2) 

n  **  2  **  ** 


(4c) 


Substituting  Equation  4  into  Equation  2  yields  a  two-term  recursion 
formula  for  updating  the  velocity: 


•n+1/2 


=  (2  -  ch/p)  ;n-l/2  _  (  2  h/p  \  n 
\2  +  ch/p /  ~  U  +  ch/p/  ~ 


(5) 


where  rn  is  the  residual  force  vector  and  is  defined  as  the  difference 
between  the  approximated  internal  and  external  force  vectors  after  n 
time  steps: 


xn  =  F.(xn)  -  F  (tn) 


(5a) 


The  displacement  is  updated  by  replacing  n  with  n+1  and  solving  for  x 
in  Equation  4a: 


n+1 


n+1  n  ,  ,  *n+l/2 

k  •  x  +  h  x 


(6) 


Substituting  Equation  5  into  Equation  6  and  using  Equation  4a  to 

•  ii"  1  /2 

eliminate  x  '  yields  a  three-term  recursion  formula  (see  Papadrakakis , 
1982)  for  updating  the  displacement  vector: 
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This  is  the  primary  formula  for  the  dynamic  relaxation  method.  It  can 
be  seen  that  the  central  difference  integration  method  results  in  an 
explicit  formulation  for  marching  through  time,  since  the  vectors  on  the 
right-hand  side  of  Equation  7  are  known.  The  vectorization  apparent  in 
Equation  7  implies  certain  data  structure  advantages  over  standard  Newton 
methods  which  are  matrix-based  formulations.  Papadrakakis  (1986)  refers 
to  this  as  a  vector  iteration  method. 

Equation  7  does  not  apply  when  n  -  0,  for  the  displacement  at  x 
is  unknown.  Therefore,  an  additional  formula  is  needed  to  start  the 
iteration  (see  Shugar,  1987)  as  follows: 


(8) 


This  is  a  consistently  derived  starting  formula  for  the  dynamic  relaxa¬ 
tion  method.  Equations  7  and  8  generally  agree  with  the  dynamic  relax¬ 
ation  formulas  given  by  Papadrakakis  (1981a)  and  Underwood  (1983).* 

It  should  be  emphasized  that  the  calculation  of  the  residual  r11  in 
Equation  7  is  computationally  straightforward  using  its  definition 
(Equation  5a).  All  that  is  required  is  to  calculate  and  sum  the  unbal¬ 
anced  forces  existing  at  each  node  point  in  the  structural  model  at  the 
known,  current  configuration  xn.  The  internal  force  calculation  for 
each  element  may  be  simply  based  on  the  conventional  linear  portion  of 
the  finite  element  stiffness  matrix,  so  long  as  the  cable  structure 
behavior  is  restricted  to  small  axial  strains.  Otherwise,  nonlinear 


*P.G.  Underwood  investigated  dynamic  relaxation  methods  for  applications 
to  nonlinear  buckling  problems  in  structural  analysis.  In  this  regard, 
the  reader  is  also  referred  to  Papadrakakis  (1981b)  and  Tong  (1986a, 
1986b). 
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stiffness  terms  may  be  derived  to  handle  large  axial  strains,  as  shown 
by  Zienkiewicz  (1977).  Large  displacements  of  cable  elements  are  per¬ 
mitted  naturally  by  the  displacement  update  formula  (Equation  7).  It  is 
not  necessary  to  use  a  geometric  stiffness  matrix  as  in  Newton-based 
methods.  The  residual  force  vector  is  calculated  on  the  updated,  local, 
or  element  level.  No  assembly  is  needed,  nor  is  it  necessary  to  compute 
the  tangent  stiffness  matrix  of  the  system  as  required  for  the 
Newton-based  procedures . 


DERIVATION  OF  OPTIMAL  ITERATION  PARAMETERS 


Since  the  artificial  inertia  and  artificial  damping  forces  appear¬ 
ing  in  the  psuedo  equation  of  motion  (Equation  2)  are  arbitrary,  one  is 
at  liberty  to  choose  any  values  for  the  parameters  p  and  c.  The  objec¬ 
tive  is  to  develop  formulas  for  these  parameters  so  that  they  may  be 
adaptively  controlled  (during  the  solution)  in  such  a  way  as  to  promote 
optimum  convergence  for  the  dynamic  relaxation  iterative  procedure. 

One  assumes  that  any  difficulty  in  achieving  convergence  is  due  to 
numerical  ill-conditioning  in  the  tangent  stiffness  matrix  K(x).  Such 
ill-conditioning  may  arise  for  valid  physical  reasons,  such  as  when 
cables  approach  slack  conditions.  Since  a  well-behaved  tangent  stiffness 
matrix  is  important  to  stability  of  the  numerical  process,  whatever  the 
desired  optimizing  formulas  are,  it  makes  sense  to  anticipate  that  they 
should  somehow  reflect  the  numerical  condition  of  the  tangent  stiffness - 
matrix  during  the  solution  process. 

To  examine  systematically  the  convergence  of  the  method,  Lynch 
(1968)  transformed  the  iterative  process  into  a  standard  eigenvalue 
problem  for  the  error  vector: 


(9) 
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Following  this  approach,  subtract  the  true  solution  x  from  Equation  7 
and  rewrite  It  using  the  error  definition  (9)  to  get  the  relationship 
between  successive  error  vectors: 


en+1  =  [B  I  -  y  B]  £n  -  a  £n_1 

(10) 

where:  B  =  2  "  +  1 

2  +  ch/p 

(10a) 

r  -  2  h2/p 

2  +  ch/p 

(10b) 

a  ■  B  -  1 

(10c) 

In  Equation  10  the  matrix  B  is  a  preconditioned  tangent  stiffness  matrix, 
B  =  W"1  K(xn)  W_T  (11) 

M  M  M  M  M 

where  the  preconditioning  matrix  V  (see  Papadrakakis ,  1986)  is: 

a.  W  =  I.  This  is  the  standard  dynamic  relaxation  scheme  which 
has  been  presented  above  in  deriving  Equation  7. 

1/2 

b.  W  =  D  .  Here,  D  is  a  diagonal  matrix  composed  of  the  entries 
on  the  main  diagonal  of  the  tangent  stiffness  matrix.  This 
method  is  also  called  diagonal  scaling  since  the  entries  on 
the  main  diagonal  of  the  coefficient  matrix,  formed  from  the 

•  1  n 

triproduct  matrix,  W  K(x  )  W  ,  are  unity. 


Letting  the  rate  at  which  the  error  vector  decays  in  the  dynamic 
relaxation  process  be  denoted  by  XpR,  then: 


n+1 

t 


(12) 


Clearly,  the  modulus  of  should  remain  less  than  unity  for  the 
dynamic  relaxation  process  to  converge  with  increasing  n. 
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Following  Lynch  (1968),  one  wishes  to  Investigate  what  Influences 
the  decay  rate  XpR.  An  equation  for  X^R  Is  obtained  by  substituting 
Equation  12  into  Equation  10.  One  can  then  obtain: 

fXDR  '  &  XDR  * 

W  *  XDR  ) 

Now,  the  standard  eigenvalue  problem  for  the  eigenpairs  (X_,en)  of  the 

D  ** 

matrix  B  is  written  as: 


n 

£  =  0  (13) 


(-  *B  l  +  £n  *  2 


(14) 


If  one  identifies  the  scalar  coefficient  of  I  in  Equation  13  with  the 
scalar  coefficient  of  I  in  Equation  14,  one  can  arrive  at: 


8  +  o 


+  XB 


0 


(15) 


In  Equation  15,  one  has  related  the  decay  rate  X^R  of  the  dynamic  relax¬ 
ation  process  to  the  eigenvalues  X_  of  the  structure's  tangent  stiffness 

D 

matrix.  This  is  the  type  of  relationship  that  was  anticipated.  Further 
more,  the  condition  of  the  structural  stiffness  matrix  can  be  tracked  by 
monitoring  its  eigenvalues  by  some  method  yet  to  be  described. 

Solving  Equation  15  for  X^R  one  obtains  the  quadratic  equation: 

XDR  -  (B  ’  *  V  XDR  +  »  -  0  <“> 


So,  for  each  of  the  N  eigenvalues  of  B  there  are  two  solutions  for  X^R: 


XDR  =  \  (8  -  y  XR)  ±  |  V(8  -  T  x„)2  -  4  a 


B  *  2 


B' 


(17) 
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There  are  three  cases  for  these  roots. 


Case  1.  The  roots  of  the  quadratic  are  complex  when: 

4a  >  (6  -  *  XB)2  (18) 


The  modulus  of  the  decay  rate  jX^j  for  this  case  is: 


The  decay  rate  modulus  is  seen  to  be  independent  of  the  eigenvalue  X_  in 

15 

the  case  of  complex  roots,  and,  therefore,  this  case  is  of  no  particular 
use. 


Case  2.  The  roots  are  real  and  equal  when: 


4  a  =  (P  -  »  Xfl)i 


(20) 


The  roots  are: 


X  =  _ l—  (2  -  £  X  ] 

DR  2  +  ch/p  V  P  Bj 


(21) 


For  this  case,  one  also  finds  that  the  square  of  the  convergence 
parameter  ch/p  has  the  following  relationship  with  Xg: 


X„  h 


(?r  ■  ^  (*  - V-j 


(22) 


Case  3.  The  roots  are  real  and  unequal  when: 


4  «  <  (3  -  r  xB)' 


(23) 
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The  larger  of  the  two  roots  Is  found  to  be: 


1 

2  +  ch/p 


From  the  definition  of  XDR  in  Equation  12,  one  can  see  that  it  is 
important  to  promote  the  smallest  value  possible  for  this  parameter,  and 
to  require  that  its  modulus  never  exceed  unity. 

With  that  in  mind,  the  behavior  of  the  modulus  |XpR|  is  studied  by 

graphing  it  versus  ch/p  using  Equation  24,  while  holding  the  parameter 

2  2 
(h  /p)  Xg  constant  (Figure  1).  One  can  see  that  for  any  value  of  (h  /p) 

X_,  the  smallest  value  of  Xn_  occurs  for  Case  2,  the  equal-roots  case. 

®  2 

Thus,  for  any  value  of  (h  /p)  X_,  the  appropriate  choice  for  ch/p  comes 

B 

from  Equation  22. 

However,  note  that  there  are  N  values  of  Xn  and,  thus,  N  values  of 

(h  /p)  Xg.  To  help  resolve  this,  the  behavior  of  |XDR|  versus  (h  /p)  XR 

is  similarly  graphed  in  Figure  2,  while  holding  constant  the  parameter 

2 

ch/p.  Not  shown  is  the  region  where  XQR  exceeds  unity  when  (h  /p)  XR  > 

4.0.  This  region  is  to  be  avoided  if  convergence  is  to  be  promoted. 

2 

Moreover,  a  strategy  for  selecting  (h  /p)  XR  is  suggested  from  this  fig¬ 
ure. 

Assume  that  for  any  B,  the  maximum  eigenvalue  lgmax  the  minimum 

eigenvalue  X_  .  have  been  calculated.  Then  a  value  of  h2/p  that  pro- 
am  in  2 

motes  convergence  can  be  calculated  by  centering  (h  /p)  XRm^n  ^ 

(h2/p)  *gmax  symmetrically  about  the  abscissa  value  of  2.  The 
appropriate  equation  is  therefore: 


4 

X_  + 
Bmax 


^Bmin 


(25) 


11 


2 

Using  this  equation  avoids  the  region  (h  /p)  >  4.0,  and  assures 

stability  of  the  dynamic  relaxation  procedure.  Furthermore,  substituting 
Equation  25  into  Equation  22  leads  to  the  optimal  choice  for  ch/p  for 
promoting  a  faster  convergence  rate  for  the  dynamic  relaxation  procedure: 


^Bmax  +  ^Bmin 


V* 


Bmax 


Bmin 


(26) 


ESTIMATING  THE  MAXIMUM  AND  MINIMUM  EIGENVALUES 

To  ensure  stability  one  needs  only  to  ensure  that  the  inequality 

|XDR|  <  1  is  maintained  during  iteration  with  the  dynamic  relaxation 

method.  Thus,  if  one  is  able  to  calculate  an  upper  bound  for  the  maxi- 

2 

mum  eigenvalue,  then  one  is  assured  that  when  its  product  with  h  /p  is 

2 

taken  to  be  less  than  four,  the  product  of  ^Bmax  with  h  /p  will  also  be 
less  than  four  and  thereby  satisfy  |XDR|  <  1.  Papadrakakis  (1981a)  sug¬ 
gests  the  following  method. 

The  Gerschgorin  theorem  is  available  to  calculate  an  upper  bound  of 
the  maximum  eigenvalue  of  a  square  matrix.  Applying  it,  one  can  write: 

N 

IW  ‘  T  l  l»ijl  (27) 

j=l 


where  the  B . .  are  the  entries  in  the  scaled  stiffness  matrix  B.  This 

ij 

states  that  the  maximum  over  all  rows  of  B  of  the  sum  of  the  entries 
(their  absolute  values)  in  each  row  is  an  upper  bound  of  ^Bmax-  An 
implementation  of  this  simple  inequality  must  involve  an  algorithm  that 
recognizes  that  the  matrix  B  is  never  formed  explicitly  in  the  dynamic 
relaxation  method.  Thus,  the  coefficients  B^  are  calculated  at  the 
element  level  to  preserve  the  advantage  of  a  vector  iteration  method. 
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Calculation  of  the  minimum  eigenvalue  is  less  critical  than  the 
calculation  of  the  maximum  eigenvalue  for  it  cannot  directly  affect  the 
stability  of  the  iteration  process.  Exact  methods  for  its  calculation 
are  not  required,  and  bounds  are  not  necessary  for  its  estimation.  The 
minimum  eigenvalue  must  only  satisfy: 

®  <  ^Bmin  <  ^Bmax  (28) 

Its  value,  along  with  that  for  the  maximum  eigenvalue,  does  affect  the 
rate  of  convergence  of  the  solution.  Numerical  experience  has  shown 
that  poor  estimates  can  adversely  affect  the  rate  of  convergence,  though 
they  cannot  directly  cause  a  solution  to  blow  up.  Good  estimates  of  the 
minimum  eigenvalue  should  therefore  be  sought. 

It  is  noted  that  Equations  27  and  28  are  only  needed  on  occasions 
when  progress  toward  a  solution  is  not  satisfactory.  Then  they  are  em¬ 
ployed  along  with  the  optimal  formulas  (Equations  25  and  26)  to  adjust 
the  solution  parameters  in  Equation  7. 


AUTOMATED  DYNAMIC  RELATION  ALGORITHM 

A  procedural  description  of  the  automated  dynamic  relaxation 
solution  method  is  given  by  the  following  algorithm: 
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KINEMATIC  DAMPING  ALGORITHM 

Cundall  (1976),  when  examining  the  application  of  explicit 
integration  methods  to  problems  in  geomechanics,  suggested  that  the 
kinetic  energy  of  the  structure  be  constantly  monitored,  and  that  when 
an  energy  peak  is  detected  all  the  current  velocities  be  set  to  zero. 
For  a  system  oscillating  in  one  mode,  the  state  of  stress  at  the  energy 
peak  would  correspond  to  the  static  equilibrium  position.  However,  for 
practical  problems  with  many  degrees  of  freedom,  the  process  must  be 
repeated  through  further  peaks,  eliminating  the  kinetic  energy  for  all 
modes,  until  the  required  degree  of  accuracy  is  obtained. 
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Using  this  method ,  the  viscous  damping  coefficient  of  Equation  2  is 
set  to  zero  and  the  equation  of  motion  becomes: 


&  +  EiC*)  -  Ee(t) 


(29) 


Integrating  Equation  29,  using  central  difference  approximations  in  the 
same  way  as  Equation  2,  we  get: 


•n+1/2  *n-l/2  .  n 

x  =  5  -  (h/p)  r 

and, 

xn+1  =  2xn  -  x-1  -  Ch2/P) 


(30) 


(31) 


Following  the  same  process  for  derivation  of  optimal  iteration  parameters 

2 

as  before,  the  optimum  value  for  (h  /p)  is: 


It  is  also  recommended  the  the  sum  of  and  ^gnax  be  replaced  by  the 

Gerschgorin  bound  in  Equation  27. 

When  the  time  increment  h  is  kept  constant  throughout  the  iteration 
process,  an  energy  peak  can  be  detected  at  each  successive  time  step, 
and  in  this  case  the  structure  is  oscillating  about  the  equilibrium  state 
with  a  constant  frequency.  To  overcome  this  difficulty,  the  time  step 
is  reduced  to  half  of  its  value  every  time  an  energy  peak  is  detected  at 
each  time  step  until  the  iteration  process  returns  to  its  normal  conver¬ 
gence.  Once  this  has  been  achieved,  the  time  increment  is  reset  to  its 
original  value.  A  flow  chart  of  the  method  is  shown  in  Figure  3. 

Papadrakakis  (1988)  believes  this  version  of  dynamic  relaxation  may 
have  better  behavior  than  the  classical  viscous  damping  version  of  dyna¬ 
mic  relaxation,  presented  here,  for  the  class  of  problems  studied.  This 
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would  seem  to  be  corroborated  by  the  success  of  Barnes  and  coworkers  (in 
the  U.K.)  who  have  used  the  kinetic  damping  version  of  dynamic  relaxation 
for  problems  involving  compliant  structures.  For  example,  flexible  risers 
were  studied  by  Soltanahmadi  and  Barnes  (1987),  prestressed  nets  and 
membranes  were  studied  by  Barnes  (1987),  and  various  other  membranes  and 
tensioned  fabric  structures  have  been  solved  by  Wakefield  (1987).  Much 
of  the  success  of  ADR  in  the  design  office  environment  is  ascribed  to 
the  simplicity  of  the  kinetic  damping  algorithm.  It  is  pointed  out  that 
the  inherent  dynamic  analogy  of  the  method  also  facilitates  understanding 
and  engenders  control  of  the  solution  procedure  by  practicing  engineers 
when  solving  complex  problems. 


RESULTS  FROM  NUMERICAL  EXPERIMENTS 

In  the  following  numerical  experiments,  results  from  the  standard 
Newton  iteration  algorithm  and  the  ADR  algorithm  applied  to  nonlinear 
static  cable  problems  are  compared.  These  results  are  abstracted  from 
Shugar  (1987).  Further  results  comparing  the  performance  of  the  ADR 
algorithm  with  existing  cable  analysis  computer  programs  are  also  pre¬ 
sented. 

A  set  of  simple,  two-dimensional  cable  test  problems  was  designed 
primarily  to  test  an  algorithm's  ability  to  seek  the  correct  static 
equilibrium  solution  from  an  initial  cable  configuration  that  is  arbi¬ 
trarily  prescribed.  The  robustness  of  the  algorithm  is  addressed  by 
deliberately  specifying  onerous  initial  nonequilibrium  configurations. 
Furthermore,  multiple  solutions  are  possible  in  some  cases  due  to  speci¬ 
fication  of  a  simple  cable  element  that  can  support  axial  compression 
force  and  has  no  bending  stiffness. 

Fixed-Span  Suspended  Cable  Problem 

The  fixed-span  suspended  cable  test  problem  consists  of  a  cable 
that  is  suspended  between  two  supports  and  is  acted  on  by  lateral,  unit 
forces  at  each  node.  The  unstrained  length  of  the  cable  is  100  wits, 
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and  the  span  is  60  units.  The  cable  length  is  uniformly  divided  into  10 
elements.  The  rigidity  of  the  cables,  EA,  is  1000.  Table  1  compares 
the  convergence  results  of  the  full  Newton  and  ADR  algorithms  for  four 
initial  configuration  cases.  The  results  indicate  that  convergence  to 
the  correct  solution  was  achieved  successfully  by  both  algorithms  for 
all  cases  except  the  third,  the  kinked  Initial  configuration  case. 

This  case  consists  of  an  initial  configuration  where  all  the 
elements  lie  on  a  straight  line  across  the  span,  and  the  two  innermost 
elements  overlap  (though  they  are  shown  parallel  with)  adjacent  ele¬ 
ments.  The  expected  equilibrium  solution  for  this  problem  is  presented 
in  Figure  4.  Also  shown  is  the  unexpected,  kinked  equilibrium  position 
found  by  the  full  Newton  algorithm  for  certain  prescribed  initial  pre¬ 
stress  forces  —  a  requirement  with  this  algorithm  to  avoid  an  initially 
singular  system.  In  this  solution  the  two  overlapped  elements  are  sus¬ 
taining  compressive  forces.  It  is  expected  that  all  algorithms  would 
find  the  expected  solution  if  provision  for  tension-only  element  behav¬ 
ior  were  made,  since  this  constraint  would  eliminate  the  alternative 
equilibrium  states.  However,  the  ADR  algorithm  requires  no  such  provi¬ 
sion,  and  furthermore,  is  independent  of  initial  prestress  force.  There¬ 
fore,  it  exhibits  greater  robustness  in  seeking  the  expected  solution. 

Cable  Snap -Through  Problem 

The  cable  snap-through  test  problem  consists  of  a  six-element 
suspended  cable  with  an  unstrained  length  of  60  units.  The  cable  is 
considered  weightless,  has  an  EA  value  of  1000,  and  is  subjected  to 
concentrated  horizontal  and  vertical  forces,  as  shown  in  Table  2.  The 
horizontally  loaded  end  support  is  free  to  move  in  the  horizontal  direc¬ 
tion,  beginning  from  either  of  two  prescribed  initial  configuration  cases, 
until  a  static  equilibrium  configuration  is  reached.  In  the  second  case, 
the  six  cable  elements  are  initially  coiled  on  top  of  one  another  such 
that  their  nodes  possess  exactly  the  same  initial  coordinate  values.  As 
the  solution  algorithm  seeks  equilibrium,  the  cable  tends  to  uncoil  and 
snap  through  from  left  to  right. 
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The  ADR  algorithm  converges  to  the  expected  equilibrium  solution, 
which  is  shown  in  Figure  5,  from  either  initial  configuration.  However, 
in  the  first  case,  the  final  calculated  span  was  short  by  9  percent. 

(This  is  an  anomalous  result;  in  the  other  13  test  problems  no  similar 
inaccuracy  occurred  with  the  ADR  algorithm.)  The  accuracy  of  the  full 
Newton  algorithm  was  good  in  this  case.  In  the  second  case,  however, 
the  full  Newton  algorithm  did  not  always  converge  to  the  expected  solu¬ 
tion;  again,  this  behavior  depended  on  the  value  of  the  required  initial 
prestress  force.  The  ADR  algorithm  is  independent  of  this  quantity,  and 
converged  to  the  expected  solution  accurately. 

TVo  unexpected,  alternative  equilibrium  states  found  by  the  full 
Newton  algorithm  are  shown  in  Figures  5b  and  5c.  Once  again  these 
states  include  elements  sustaining  compression  forces. 

Varying-Span  Suspended  Cable  Problem 

The  varying-span  suspended  cable  test  problem  consists  of  a  hori¬ 
zontally  suspended  cable  with  one  of  the  two  supports  free  to  slide 
horizontally.  The  cable  is  acted  upon  by  uniform  lateral  load  forces 
and  a  concentrated,  horizontal  reaction  force,  both  of  which  cause  one 
support  to  slide  until  the  system  reaches  equilibrium.  The  unstrained 
cable  is  200  feet  long  and  is  uniformly  subdivided  into  10  elements. 

This  test  problem  was  used  by  Webster  (1979)  in  a  similar  study  of 
solution  algorithms  for  ocean  cable  systems. 

The  three  initial  configuration  cases  studied  are  shown  in  Table  3 
along  with  the  problem  parameters.  In  the  first  two  cases,  the  cable  is 
initially  aligned  horizontally  along  the  span.  Case  1  is  labeled  taut 
because  its  initial  configuration  is  such  that  the  cable  will  remain  in 
tension  throughout  the  solution  process.  Case  2  is  labeled  slack  because 
the  cable  must  pass  from  a  compression  state  to  a  tension  state  during 
the  solution,  thus  traversing  a  zero- tens  ion  state.  In  Case  3,  a  com¬ 
pletely  sagged  cable  is  represented  wherein  all  elements  initially  lie 
along  a  common  vertical  line;  the  cable  should  otherwise  remain  in  ten¬ 
sion  throughout  the  solution  process. 
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The  ADR  algorithm  converged  to  the  expected  equilibrium  solution, 
shown  in  Figure  6,  in  all  three  cases.  However,  in  Case  2,  the  full 
Newton  algorithm  did  not  always  converge  to  this  solution;  its  behavior 
was  like  that  already  discussed  for  the  previous  problems.  Whether  the 
expected  solution  was  found  or  not  depended  on  the  value  of  the  prescribed 
prestress  force.  Webster  (1979)  showed  that  the  full  Newton  algorithm 
failed  to  converge  for  this  test  case.  Here,  it  was  learned  that  it  may 
or  may  not  converge  to  the  expected  solution.  In  general,  however,  these 
results  are  in  agreement  with  those  of  Webster's  concerning  the  unreli¬ 
ability  of  the  full  Newton  algorithm  for  this  class  of  problems. 

The  ADR  algorithm  is  more  robust  than  the  full  Newton  algorithm  in 
seeking  the  expected  equilibrium  solution  in  the  varying-span  suspended 
cable  problem.  Once  again,  this  behavior  is  correlated  with  whether  or 
not  a  slack  cable  configuration  arises  in  the  solution  process. 

Also  the  results  from  Cases  1  and  3  show  that  the  ability  of  either 
algorithm  to  converge  is  not  dependent  on  the  degree  of  cable  sag  present 
in  the  initial  configuration.  Thus,  for  numerical  conditioning,  it  seems 
to  matter  whether  cable  slackness  (zero  tension)  is  present,  but  not  to 
what  degree  cable  sag  is  present. 

Comparison  with  SEADYN  and  SEASTAR 

In  Table  4,  the  results  of  the  preceding  numerical  experiments  are 
summarized  for  comparison  with  the  results  of  further  experiments  on  the 
same  set  of  test  problems  using  two  different  nonlinear,  displacement- 
based  finite  element  programs  which  are  currently  in  practice  for  static 
cable  analysis  applications,  SEADYN  and  SEASTAR.  Both  of  these  programs 
were  evaluated  over  the  set  of  14  (problem  7  is  actually  four  similar 
problems)  test  problems,  but  did  not  perform  quite  as  well  as  the  ADR 
solution  algorithm. 
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The  results  from  SEASTAR,  using  its  Newton-Raphson  solution  option* 
(the  default  option),  are  very  similar  to  the  results  presented  earlier 
for  the  test  code  written  for  this  study  which  uses  the  same  solution 
algorithm  (thus  this  similarity  is  to  be  expected).  SEASTAR  was  able  to 
converge  to  the  expected  stable  equilibrium  solution  for  each  problem. 
However,  it  was  also  able  to  converge  to  unexpected  cable  equilibrium 
positions  as  well.  The  difference  is  a  consequence  of  having  specified 
different  values  of  certain  input  parameters.  The  parameters  with  the 
most  noticeable  influence  were  the  initial  pretension  of  the  cable 
elements,  the  load  step  or  increment  size,  and  whether  the  cable 
elements  were  allowed  to  sustain  compression. 

When  a  load  step  size  of  1.0  was  prescribed  (i.e.,  all  the  load  is 
applied  at  once),  SEASTAR  converged  to  the  expected  equilibrium  solution. 
However,  in  test  problem  3,  application  of  the  external  load  in  increments 
instead  of  all  at  once,  a  common  and  reasonable  strategy  for  nonlinear 
problems,  resulted  in  convergence  to  only  metastable  solutions.  Figures 
7a  and  7b  show  this  outcome  when  the  total  load  was  applied  in  two  and 
four  increments,  respectively.  The  unexpected  solution  obtained  by 
SEASTAR  for  test  problem  6  was  similar  to  the  profile  obtained  before 
(see  Figures  4  and  5)  and  is  omitted. 

In  test  problem  10,  application  of  the  total  load  in  one  step  resulted 
in  convergence  to  the  expected  equilibrium  solution,  as  shown  in  Figure 
8a.  Unfortunately  this  behavior  seems  unpredictable,  for  application  of 
the  load  in  four  increments  and  20  increments  again  resulted  in  convergence 
to  only  metastable  solutions,  as  shown  in  Figures  8b  and  8c,  respectively. 
Prescribing  a  zero  pretension  load  is  not  necessarily  the  answer  in  this 
test  problem  as  indicated  by  the  result  in  Figure  8d. 


♦SEASTAR  is  a  modified  version  of  a  computer  program  called  ANSR  III 
which  was  written  for  general  purpose,  nonlinear  finite  element, 
structural  analysis  applications.  As  such,  it  has  additional  advanced 
options  for  solving  nonlinear  static  problems  (Simons  and  Powell,  1982). 
These  options  must  be  used  with  caution  and  require  experience  for  they 
are  based  on  intricate  modifications  of  the  Newton-Raphson  algorithm 
aimed  at  overcoming  its  inherent  unreliability  in  structural  buckling 
problems. 
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The  evolution  of  the  unexpected  solution  (Figure  8d)  is  presented 
in  Figure  8e,  which  shows  four  intermediate  plots  depicting  initial 
progress  toward  the  converged  unexpected  solution.  Load  increments  1 
through  5  applied  5  percent  of  the  total  load  of  each  increment.  From  5 
percent  to  100  percent  of  the  load  at  increments  of  5  percent,  the  unex¬ 
pected  solution  was  locked  in.  Thus,  the  unexpected  solution  locked  in 
early. 

The  results  from  SEADYN's  viscous  relaxation  solution  algorithm 
(called  VRR)  are  perfect  except  for  test  problem  6  where  only  the 
unexpected,  or  metastable  equilibrium  solution  was  obtained. 


SUMMARY  AND  CONCLUSIONS 

The  automated  dynamic  relaxation  (ADR)  algorithm  possesses  some 
attractive  theoretical  features  for  solving  ill-conditioned  numerical 
models  of  static  ocean  cable  structures.  A  major  feature  is  the  ability 
to  control  the  stability  of  the  solution  process  automatically.  Typi¬ 
cally,  when  seeking  a  prestressed,  equilibrium  configuration  from  an 
arbitrary  starting  configuration,  the  structure  stiffness  matrix  will 
exhibit  pathological  behavior  when  slack  conditions  in  the  cable  struc¬ 
ture  arise  during  the  solution  process.  The  alternative  Newton  algorithm 
has  no  provision  for  controlling  the  condition  of  the  structure  stiffness 
matrix  and  may  therefore  be  defeated  by  slack  cable  behavior. 

The  Newton  algorithm  requires  that  a  cable  pretension  force  be 
arbitrarily  prescribed  to  avoid  a  singularity  condition  in  the  cable 
structure’s  stiffness  matrix  at  the  beginning  of  the  solution  process. 
This  study  has  shown  that  robustness  is  adversely  affected  by  this  re¬ 
quirement.  Convergence  to  the  expected  equilibrium  solution  sometimes 
depended  on  the  value  of  the  arbitrarily  prescribed  pretension  force 
when  multiple  equilibrium  solutions  were  present.  Input  data  require¬ 
ments  that  affect  convergence  in  this  way  run  contrary  to  the  goal  of  a 
foolproof  solution  algorithm. 
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The  ADR  algorithm  proved  to  be  more  robust  than  the  full  Newton 
algorithm  and  the  solution  algorithms  in  SEASTAR  and  SEADYN  for  the 
class  of  problems  studied.  It  converged  to  the  expected  static  equili¬ 
brium  configuration  in  every  test  problem.  Conversely,  solution  algo¬ 
rithms  based  on  the  full  Newton  method  sometimes  failed  to  converge  to 
the  expected  solution,  and  instead  converged  to  unexpected  and  sometimes 
meaningless  alternative  equilibrium  states  in  those  cases  where  the  solu¬ 
tion  process  encountered  slack  cable  behavior. 


ACKNOWLEDGMENT 

This  report  was  partially  prepared  while  the  author  was  visiting 
the  Ocean  Engineering  Laboratory,  Department  of  Mechanical  Engineering, 
University  of  California,  Santa  Barbara,  under  the  auspices  of  the  Office 
of  Naval  Research.  Jim  Cox,  Steve  Karnoski,  and  Warren  Bartel,  all  of 
NCEL,  provided  excellent  support  in  the  preparation  of  the  numerical 
studies . 


REFERENCES 

Barnes,  M.R.  (1987).  "Form-finding  and  analysis  of  prestressed  nets  and 
membranes,"  Non-Conventional  Structures,  vol  1,  B.H.V.  Topping,  Ed., 
Civil-Comp  Press,  Edinburgh,  Dec  1987,  pp  327-338. 

Cox,  J.V.  (1982).  STATMOOR  -  a  single-point  mooring  static  analysis 
program,  Naval  Civil  Engineering  Laboratory,  Technical  Note:  N-1634. 

Port  Hueneme,  CA,  Jun  1982. 

Crisfield,  M.A.  (1986).  Finite  elements  and  solution  procedures  for 
structural  analysis,  Chapter  7.  Swansea,  U.K.,  Pineridge  Press,  1986. 

Cundall,  P.A.  (1976).  "Explicit  finite  difference  methods  in  geomechanics," 
Conference  on  Numerical  Methods  in  Geomechanic,  Blacksburg,  VA,  Jun  1976. 


22 


Kamoski,  S.,  and  Palo,  P.  (1986).  FLEETMOR  validation  report,  Naval 
Civil  Engineering  Laboratory,  Memorandum  to  files.  Port  Hueneme,  CA, 

Oct  1986. 

Liu,  F.C.  (1977).  Establishment  of  initial  configurations  in  lumped 
parameter  simulations  of  underwater  cable  dynamics.  Naval  Civil 
Engineering  Laboratory,  Memorandum  to  files.  Port  Hueneme,  CA, 

May  1977. 

Lynch,  R.D.  (1968).  The  application  of  dynamic  relaxation  to  the  finite 
element  method  of  structural  analysis,  Ph.D.  Dissertation,  University  of 
Notre  Dame,  IN,  Jul  1968. 

Naval  Facilities  Engineering  Command  (1985).  Fleet  moorings.  Design 
Manual  DM-26.5.  Alexandria,  VA,  Jun  1985. 

O'Brien,  T. ,  and  Francis,  A.J.  (1964).  "Cable  movements  under  two- 
dimensional  loads,"  Journal,  Structural  Division,  American  Society  of 
Civil  Engineers,  vol  19,  no.  ST3,  Jun  1964,  pp  89-123. 

Palo,  P.,  and  Kamoski,  S.  (1986).  FLEETMOR  -  Static  mooring  analysis 
model  users  manual.  Naval  Civil  Engineering  Laboratory,  Memorandum  to 
files.  Port  Hueneme,  CA,  Jul  1986. 

Papadrakakis ,  M.  (1981a).  "A  method  for  the  automatic  evaluation  of  the 
dynamic  relaxation  parameters,"  Computer  Methods  in  Applied  Mechanics 
and  Engineering,  1981,  pp  35-48. 

_  (1981b).  "Post-buckling  analysis  of  spatial  structures  by 

vector  iteration  methods,"  Computers  and  Structures,  14,  1981, 
pp  393-402. 


_  (1982).  "A  family  of  methods  with  three-term  recursion 

formulae,"  International  Journal  for  Numerical  Methods  in  Engineering, 
vol  18,  1982,  pp  1785-1799. 


23 


_  (1986).  "Accelerating  vector  interation  methods,"  Journal  of 

Applied  Mechanics,  vol  53,  Jun  1986,  pp  291-297. 

_  (1988).  Personal  communication,  Department  of  Civil  Engineering, 

National  Technical  University,  Athens,  Greece,  10  Mar  1988. 

Pawsey,  S.J.,  and  Nour-Omid,  S.  (1988).  Personal  communication  and 
SEASTAR  users  manual,  PMB  Engineering,  Inc.,  San  Francisco,  CA,  Feb  1988. 

Peyrot,  A.H.,  and  Goulois,  A.M.  (1979).  "Analysis  of  cable  structures," 
Computers  and  Structures,  vol  10,  1979,  pp  805-813. 

Peyrot,  A.H.  (1980).  "Marine  cable  structures,"  Journal,  Structural 
Division,  American  Society  of  Civil  Engineers,  vol  106,  no.  ST12,  Dec  1980. 

_  (1989).  "Cables  and  flexible  pipes  in  contact  with  the  sea 

floor,"  Steel  Structures,  in  Proceedings  of  Structures  Congress  *89, 

S.B.  Iffland,  Ed.,  New  York,  N.Y.,  American  Society  of  Civil  Engineers, 

1989,  pp  31-40. 

Shields,  D.R.,  and  Zueck,  R.F.  (1984).  "Deepwater  semisubmersible  motion 
simulation,"  paper  presented  at  OCEANS  '84,  Washington,  DC,  1984. 

Shugar,  T.A.  (1987).  A  numerical  study  of  automated  dynamic  relaxation 
for  nonlinear  static  tensioned  structures,  Naval  Civil  Engineering 
Laboratory,  Technical  Note  N-1778.  Port  Hueneme,  CA,  Oct  1987. 

Shugar,  T.A.,  and  Armand,  J.-L.  (1987).  Computational  methods  in  ocean 
structural  engineering  -  a  review,  Naval  Civil  Engineering  Laboratory, 
Technical  Note  N-1768.  Port  Hueneme,  CA,  May  1987. 

Shugar,  T.A.  (1988).  "Automated  dynamic  relaxation  solution  method  for 
ocean  cable  structures,"  in  Proceedings,  Ocean  Structural  Dynamics 
Symposium  '88,  Oregon  State  University,  Corvallis,  OR,  13-14  Sep  1988, 
pp  186-205. 


24 


Simons,  J.W. ,  and  Powell,  6.H.  (1982).  Solution  strategies  for 
statically  loaded  nonlinear  structures,  Earthquake  Engineering  Research 
Center,  UCB/EERC -82/22,  Berkeley,  CA,  Nov  1982. 

Soltanahmadi,  A.,  and  Barnes,  M.R.  (1987).  "Dynamic  response  and 
modeling  of  flexible  risers,"  Non-Convent ional  Structures,  vol  1, 

B.H.V.  Topping,  Ed.,  Civil-Comp  Press,  Edinburgh,  Dec  1987,  pp  97-105. 

Tong,  P.  (1986a).  "An  adaptive  dynamic  relaxation  method  for  static 
problems,"  International  Journal  for  Computational  Mechanics,  1,  1986, 
pp  127-140. 


_ .  (1986b).  "Nonlinear  problems  by  adaptive  dynamic 

relaxation  methods,"  in  Proceedings  of  the  Army  Symposium  on  Solid 
Mechanics,  MTL-MS-86,  and  U.S. -China  Workshop  on  FEM,  Chandu-China, 

1986. 

Underwood,  P.G.  (1983).  "Dynamic  relaxation  -  a  review,"  Computational 
Methods  for  Transient  Dynamic  Analysis,  T.  Belytshcki  and  T.J.R.  Hughes, 
editors.  North  Holland,  Amsterdam,  1983,  Chapter  5. 

Wakefield,  D.S.  (1987).  "Practical  numerical  modeling  of  complex 
structures,"  Non-Conventional  Structures,  vol  2,  B.H.V.  Topping,  Ed., 
Civil-Comp  Press,  Edinburgh,  Dec  1987,  pp  75-80. 

Webster,  R.L.  (1976).  An  application  of  the  finite  element  method  to 
the  determination  of  nonlinear  static  and  dynamic  responses  of  under 
water  cable  structures,  Ph.D.  thesis,  Cornell  University,  Ithaca,  NY, 
1976. 


_  (1977).  "Finite  element  analysis  of  deep  sea  moors  and  cable 

systems,"  American  Society  of  Civil  Engineers  Fall  Convention,  Oct  1977. 
(Preprint  3033) 


25 


_  (1979).  "On  the  static  analysis  of  structures  with  strong 

geometric  nonlinearity,"  paper  presented  at  the  Second  International 
Conference  on  Computational  Methods  in  Nonlinear  Mechanics,  University 
of  Texas,  Austin,  TX,  26-29  Mar  1979. 

Webster,  R.L. ,  and  Palo,  P.A.  (1982).  SEADYN  user's  manual.  Naval  Civil 
Engineering  Laboratory,  Technical  Note  N-1630.  Port  Hueneme,  CA,  Apr  1982. 

Webster,  R.L.  (1984).  Personal  communication,  Consulting  Engineer,  Brigham 
City,  UT,  5  Jan  1984. 

Zienkiewicz,  O.C.  (1977).  The  finite  element  method.  Chapter  19,  3rd 
Ed.,  London,  England,  McGraw-Hill  Book  Co.,  1977. 


26 


Table  1.  Numerical  Results  of  Fixed-Span  Suspended  Cable 


Case 

No. 


Initial 

Configuration 


Algorithms 


Full 

Newton 


ADR 


1  Rectangle 


Converged 


Converged 


Converged 


Converged 


(all  elements  are  colinear) 
4  Saw  tooth 


Converged  Converged 
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Table  2.  Numerical  Results  -  Cable  Snap  Through 


Case 

No. 


Initial 

Configuration 


Algorithms 


Full 

Newton 


ADR 
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Table  3.  Numerical  Results  -  Varying-Span  Suspended  Cable 

Case  Initial 

No.  Configuration 

1  Taut 


H 

Converged  Converged 


ilium 


2  Slack 


Can  converge 
to  unexpected 


-O—O—O— O—O— O—O— O—O  to  unexp 

|  1  j  j  j  j  j  \/r>r 


w 


Converged 


Sag 


(all  elements 
are  colinear) 


L  =  200  ft 
H  =  5.77  lb 
w  =  0. 1  lb/ft 
EA  =  1  x  10*  lb 


Converged 


Converged 


Table  4.  Summary  of  Numerical  Convergence  Performance 
for  Various  Algorithms 


Test 

Problem 

Initial 

Configuration 

Newton- 

Raphson 

ADR 

SEADYN 

VRR 

SEASTAR 

N-R 

1 

Rectangle,  suspended  cable 

E 

E 

E 

E 

2 

Triangle,  suspended  cable 

E 

E 

E 

E 

3 

Kink,  suspended  cable 

E,U 

E 

U 

E,U 

4 

Sawtooth,  suspended  cable 

E 

E 

E 

E 

5 

Triangle,  cable  snap-through 

E 

E 

E 

E 

6 

Coil,  cable  snap-through 

E,U 

E** 

E 

E,U 

7 

Inverted  sag,  four  mooring  cables 

E 

E 

E 

E 

8 

Taut  mooring  cable 

E 

E 

E 

E 

9 

Zero  sag,  varying  cable  span 

E 

E 

E 

E 

10 

Slack,  varying  cable  span 

E,U  . 

E 

E 

E,U 

11 

Max  sag,  varying  cable  span 

E 

E 

E 

E 

Key:  U  =  Converged  to  an  unexpected  or  metastable  equilibrium  solution. 
E  =  Converged  to  the  expected,  stable  equilibrium  solution. 

**  =  Converged  but  with  9%  error  in  final  span--anomalous  result. 
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Figure  1.  Dynamic  relaxation  --  relation  between 
convergence  and  ch/p. 
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Figure  3.  Row  chart  for  dynamic  relaxation  with  kinetic  damping  (Papadrakakis,  1988). 
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Figure  4.  Expected  and  unexpected  equilibrium  configuration  for 
initially  kinked  cable. 
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Figure  5.  Expected  and  unexpected  equilibrium  configurations  for 
initially  coiled  cable. 
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Initial  load 
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Figure  8d.  Varying-span  cable,  case  2  -  slack.  Initial  load 
load  step  =  0.05. 
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Figure  8e.  Varying-span  cable,  case  2  -  slack.  Initial  load 
load  step  =  0.05. 


Appendix 


DISCUSSION  OF  EXISTING  NUMERICAL  PROCEDURES 
FOR  MOORING  AND  OCEAN  CABLE  SYSTEMS 


The  design  of  mooring  systems  is  generally  governed  by  static  load 
behavior.  The  mooring  problem  described  here  is  a  static  force  problem 
in  the  horizontal  plane  for  which  the  conditions  of  static  equilibrium 
are  imposed  by  three  equations  of  equilibrium.  Single-point  mooring 
problems  are  statically  determinant  and  multipoint  mooring  problems  are 
statically  indeterminant.  For  the  single-point  mooring  problem,  there 
are  available  the  three  equations  of  equilibrium  in  the  plane  to  solve 
for  the  three  unknowns;  the  ship  heading  angle,  the  hawser  force,  and 
the  hawser  heading  angle.  This  problem  can  therefore  be  solved  directly 
without  resorting  to  the  force-displacement  relation  and  compatibility 
condition  for  the  compliant  mooring  line.  For  multipoint  mooring  prob¬ 
lems,  each  additional  mooring  leg,  beyond  the  first,  contributes  two 
more  unknowns  (either  two  rectangular  components  of  the  restoring  force 
or  the  restoring  force  and  its  angle)  to  the  mooring  problem.  For  an 
m-point  mooring  problem,  the  order  of  redundancy  is  therefore  2(m-l). 
Thus,  2(m-l)  compatibility  equations  are  required,  in  addition  to  the 
three  equations  of  equilibrium,  to  solve  the  statically  indeterminant 
problem. 

The  required  compatibility  equations  are  provided  by  the  displace¬ 
ment  constraints  imposed  by  the  rigid  ship  deck  and  by  the  compliant 
mooring  lines.  In  numerical  procedures,  these  constraints  are  formed 
from  two  sets  of  input  data  for  each  mooring  leg;  the  coordinates  of  the 
anchor  point  and  chock,  where  the  mooring  line  is  assumed  to  terminate, 
and  the  force-deflection  relationship  of  the  mooring  legs.  Thus,  the 
relative  distances  among  the  chocks  are  known  to  the  computer  program 
and  remain  fixed  in  the  horizontal  plane  according  to  the  first  set  of 
data,  and  the  displacement  of  the  chocks  relative  to  the  anchor  points 
can  be  calculated  using  the  second  set  of  data.  From  this  information, 
the  displacement  of  the  ship  and  the  displacements  of  the  mooring  lines 
are  calculated  such  that  compatibility  is  satisfied.  The  mooring  line 
forces  which  act  on  the  ship  are  calculated  consistent  with  the  ship 
displacement. 

The  Naval  Facilities  Engineering  Command  (NAVFAC)  promogates  a 
catenary  mooring  analysis  computer  program  called  MOORING  to  calculate 
the  position  of  a  fleet  moored  ship  and  the  forces  in  the  mooring  lines 
due  to  wind  and  current  loads.  The  analysis  and  design  procedure,  as 
well  as  this  program,  are  described  in  NAVFAC  Design  Manual  DM-26.5  (see 
the  list  of  references).  This  program  Is  a  simple  microcomputer  program 
for  both  fleet  and  fixed  mooring  design.  It  is  written  in  Microsoft 
GBASIC  language  and  operates  under  the  CPM  operating  system.  It  is 
currently  being  rewritten  in  FORTRAN.  There  are  several  attractive 
user-oriented  features  in  the  program,  and  it  is  apparently  convenient 
for  the  design  of  standard  moorings. 


A-l 


MOORING  pertains  to  a  specific  anchor  leg  definition  which  is  lim¬ 
ited  to  a  certain  configuration  and  to  two  different  types  of  materials. 

An  anchor  leg  configuration  may  include  a  hawser  from  the  ship  to  a  sur¬ 
face  buoy,  a  riser  chain  from  the  buoy  to  a  ring  or  sinker,  and  up  to 
three  ground  leg  chains  extending  from  the  ring  to  three  anchor  points. 

Any  combination  or  portion  of  this  pattern  may  be  prescribed  for  each 
leg  of  the  mooring  system.  Length,  submerged  weight,  and  breaking 
strength  of  each  line  component  is  prescribed,  and  zero  is  entered  where 
any  generic  component  in  the  anchor  leg  definition  is  to  be  omitted. 

The  nonlinear  load-displacement  curve  for  each  prescribed  anchor  chain 
configuration  is  precomputed  according  to  catenary  equations  and  saved 
for  subsequent  use  in  solving  the  nonlinear  static  equilibrium  equations 
governing  the  offset  position  of  the  moored  vessel.  Both  the  mooring 
leg  configuration  prescriptions  and  the  catenary  solutions  are  accom¬ 
plished  by  running  a  subprogram  CATZ. 

The  solution  to  the  nonlinear  static  equilibrium  equations  provides 
the  position  of  the  moored  vessel  in  terms  of  three  degrees  of  freedom; 
they  are  the  surge,  sway,  and  yaw  displacements  relative  to  an  assumed, 
prescribed  initial  position  for  the  vessel.  The  Newton-Raphson  solution 
method  is  employed.  It  has  been  noted  that  on  occasion,  convergence  to 
the  solution  does  not  occur,  and  the  computation  is  automatically  termi¬ 
nated  after  a  preset  number  of  iterations.  It  is  interesting  that, 
according  to  the  documentation  in  DM-26.5,  this  happens  when  the  mooring 
system  is  exceptionally  slack. 

Calculations  with  MOORING  were  compared  to  calculations  with 
another  ship  mooring  program,  FLEETMOR  (Palo  and  Karnoski,  1986),  and 
both  were  compared  to  measurements  of  moored  ships  made  in  the  Carquinez 
Straits  just  northeast  of  San  Francisco  (Karnoski  and  Palo,  1986).  The 
results  showed  that  FLEETMOR  was  more  accurate  than  MOORING.  One  reason 
is  that  FLEETHOR  allows  the  wind  and  current  loads  on  the  moored  vessel 
to  vary  according  to  changing  position  of  the  vessel  during  the  solution 
process,  whereas  MOORING  assumes  the  loads  to  remain  constant  and  inde¬ 
pendent  of  vessel  displacement. 

Like  MOORING,  FLEETMOR  is  also  intended  to  be  an  easy  to  use  micro¬ 
computer  program  for  analysis  of  standard  ship  mooring  systems.  It  com¬ 
putes  the  static  offset  position  of  the  moored  ship  due  to  wind,  current, 
and  wave  drift  force  loads.  It  models  up  to  24  mooring  legs  with  each 
leg  comprised  of  up  to  three  cable  segments  and  a  buoy  or  sinker.  Line 
stretch  is  included  to  simulate  the  elongation  of  synthetic  lines.  Either 
surface  or  subsurface  buoys  may  be  simulated. 

As  in  MOORING,  the  force-deflection  relationship  for  each  prescribed 
mooring  leg  subsystem  in  FLEETMOR  is  precomputed  using  two-dimensional 
catenary  formulas  and  stored  on  disc  file  for  later  retrieval  by  the 
program's  solution  routine.  Mooring  leg  configurations  anticipate  and 
are  limited  to  standard  Navy  ship  mooring  practice.  Arbitrary  mooring 
configurations  for  buoys,  ships,  platforms,  risers,  etc.,  are  not  pro¬ 
vided  for  in  these  programs.  No  cable  drag  forces  are  allowed,  and 
therefore  so-called  mooring  dominated  systems  cannot  be  analyzed  with 
these  programs.  MOORING  and  FLEETMOR  are  not  intended  for  general  pur¬ 
pose  mooring  system  analysis. 
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FLEETMOR  is  coded  in  FORTRAN  using  the  Professional  FORTRAN  com¬ 
piler  for  the  IBM-AT  (with  math  coprocessor)  version  and  the  Supersoft 
FORTRAN  compiler  for  the  Zenith  version.  It  also  has  graphics  routines 
written  in  BASIC.  Program  modification  must  be  anticipated  if  other 
FORTRAN  compilers  are  used  when  porting  the  program. 

Like  MOORING,  FLEETMOR  includes  user-oriented  editing  features. 
However,  unlike  MOORING,  FLEETMOR  also  computes  the  applied  wind,  cur¬ 
rent,  and  wave  drift  force  loads  from  user-prescribed  environmental 
data,  and  for  user-prescribed  Navy  vessels.  Accomodating  the  fact  that 
these  forces  vary  with  the  ship's  heading  relative  to  the  direction  of 
the  environment  adds  substantially  to  the  degree  of  nonlinearity  in  the 
static  equations  of  equilibrium  to  be  solved.  As  such,  FLEETMOR  solves 
a  more  difficult  problem  than  MOORING,  but  it  yields  a  more  accurate 
prediction  for  the  static  offset  position.  That  is,  to  the  extent  that 
the  ship's  assumed  initial  heading  is  different  from  its  final  heading 
in  the  static  offset  position,  FLEETMOR' s  predictions  of  the  offset 
position  and  mooring  line  tensions  will  be  more  accurate  than  MOORING's 
predictions.  In  this  sense,  a  solution  obtained  by  MOORING  can  be  re¬ 
garded  as  only  an  approximation  to  a  solution  obtained  by  FLEETMOR  given 
the  same  assumed  initial  headings.  The  differences  should  be  small  for 
standard  ship  mooring  applications  which  tend  to  be  relatively  stiff  or 
taut  multipoint  systems,  but  for  more  flexible  systems  or  for  single¬ 
point  mooring  systems,  the  differences  could  be  significant. 

The  solution  method  for  the  nonlinear  static  equations  in  FLEETMOR 
is  based  on  a  simple  bisection  technique  (method  of  false  position).  It 
iteratively  finds  the  heading  angle  of  the  ship  that  corresponds  to  a 
zero  value  of  the  residual  moment  acting  on  the  ship.  This  method  has  a 
slower  rate  of  convergence  than  the  Newton -Raphs on  solution  method  used 
in  MOORING.  However,  it  may  also  be  more  robust  for  difficult  mooring 
problems . 

The  solution  routine  in  FLEETMOR  is  in  fact  derived  from  the 
bisection  solution  routine  in  another  microcomputer  program  called 
STATMOOR  (Cox,  1982).  This  program  was  written  for  single-point  mooring 
analysis  of  Navy  ships.  It  has  no  provision  for  mooring  leg  data  and 
does  not  compute  a  static  offset  position  for  the  moored  ship.  Thus, 
STATMOOR  is  not  intended  to  be  a  mooring  analysis  program  per  se. 

In  single-point  moorings  the  equilibrium  heading  of  the  ship 
depends  only  on  the  directions  of  the  environmental  loads  and  is  not  a 
function  of  the  mooring  leg  restoring  force.  STATMOOR  calculates  this 
heading  in  the  same  way  that  FLEETMOR  does  for  multipoint  mooring  cases, 
the  only  difference  being  that  no  restoring  force  from  the  mooring  leg 
exists  in  the  residual  moment  equation.  The  equilibrium  heading  found 
by  the  bisection  routine  corresponds  to  equilibrium  among  the  forces 
acting  on  the  ship  due  to  environmental  loading,  which  are  also  coded  to 
vary  as  a  function  of  heading,  as  in  FLEETMOR. 

Once  the  equilibrium  heading  is  calculated,  STATMOOR  calculates  the 
hawser  force  and  heading  necessary  to  hold  the  static  offset  position, 
but  not  the  coordinates  of  the  position  itself.  FLEETMOR,  has  the 
necessary  mooring  leg  data  to  also  calculate  the  static  offset  position. 
In  fact,  FLEETMOR  contains  the  STATMOOR  code  and  defers  to  this  section 
of  the  program  when  a  user  flags  the  single-point  mooring  option.  It 
then  computes  the  equilibrium  heading  and  hawser  force  as  well  as  the 
static  offset  position. 
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For  deep  water  mooring  applications  where  the  drag  force  on  the 
mooring  line  becomes  a  consideration  in  the  equilibrium  of  the  system, 
and  for  general-purpose  mooring  and  ocean  cable  analysis,  none  of  the 
above  Navy  ship  mooring  computer  programs  are  suitable.  A  general- 
purpose  mooring  analysis  capability  is  required  for  p  actical  solutions 
in  these  types  of  problems.  Two  general-purpose  computer  programs  are 
in  use  for  Navy  applications,  SEADYN  (Webster,  1976,  and  Webster  and 
Palo,  1982)  and  SEASTAR  (Pawsey  and  Nour-Omid,  1988).  These  computer 
programs  are  displacement -based  finite  element  programs,  and  it  is  this 
technology  that  provides  for  general-purpose  capability.  Cables  are 
modeled  as  a  series  of  one-dimensional  finite  elements  which  may  be  sub¬ 
jected  to  hydrostatic  or  hydrodynamic  loads,  and  they  may  be  assembled 
into  an  entire  mooring  leg  or  arbitrarily  into  many  legs,  and  they  may 
be  networked  in  arbitrary  geometrical  configurations.  SEADYN  and  SEASTAR 
would  seem  more  suitable  for  mooring  and  ocean  cable  system  analysis 
than  alternative  commercially  available  general-purpose,  nonlinear  finite 
element  programs  (MARC,  ADINA,  ABAQUS,  ANSYS,  MSCNASTRAN)*  because  they 
have  been  specified  to  ocean  engineering  problems.  Their  nonlinear  analy¬ 
sis  features  are,  however,  too  extensive  to  be  appropriate  for  current 
microcomputer  implementation. 

One  computer  program  called  SOUPLE  (Peyrot,  1980)  is  implemented  on 
IBM-AT  microcomputers  and  provides  general-purpose  mooring  and  cable 
analysis  capability.  This  program  is  based  on  the  catenary  formulas 
given  by  O'Brien  and  Francis  *.1964)  for  use  in  modeling  mooring  legs  as 
in  MOORING  and  FLEETMOR,  but  the  catenary  force-deflection  relationship 
is  cast  (Peyrot  and  Goulois,  1979)  as  a  flexibility  relationship  as 
described  in  the  theory  of  matrix  analysis  of  structures.  This  rela¬ 
tionship  is  then  simply  inverted  to  obtain  a  stiffness  relationship  so 
that  it  behaves  like  a  super  cable  finite  element.  This  is  so  that  the 
standard,  direct  stiffness  algorithm  may  be  used  to  assemble  arbitrary 
networks  of  these  elements,  as  well  as  other  structural  elements,  into  a 
structural  model  of  a  mooring  or  cable  system.  Depending  on  the  problem, 
either  a  single  catenary  element  or  a  series  of  catenary  elements  are 
used  to  model  each  mooring  leg.  Thus,  SOUPLE  is  a  displacement-based 
structural  analysis  computer  program  which  has  been  specified  to  ocean 
engineering  problems  (Peyrot,  1980  and  1989). 

Hydrodynamic  loads  are  distributed  to  the  cable  elements  in  such  a 
way  that  the  catenary  equations  (classically  associated  with  gravity 
loading  only)  upon  which  they  are  based  remain  valid.  In  a  sense,  the 
loads  are  treated  as  an  artificial  gravity  force  component  distributed 
along  a  cable,  but  in  the  direction  of  the  normal  drag  force  component 
instead  of  in  the  vertical  direction  typical  of  conventional  catenary 
mooring  legs.  The  tangential  drag  force  component  is  neglected. 


*Some  of  these  programs  (ABAQUS,  for  example)  do  provide  for  ocean 
cable  analysis  including  hydrodynamic  loads  and  spectral  analysis 
capability.  They  should  be  considered  where  they  are  available 
when  solving  ocean  structural  analysis  problems  because  they  are 
user -oriented,  commercially-tested  products. 
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SOUPLE  does  not  provide  for  variable  hydrodynamic  forces  on  a 
moored  ship  as  does  FLEETMOR,  and  therefore  it  would  not  be  as  appro¬ 
priate  in  ship  mooring  problems  where  the  specification  of  such  forces 
is  important.  However,  it  can  be  recommended  for  use  in  other  or  more 
advanced  problems  when  the  larger  general-purpose  programs  SEADYN  and 
SEASTAR  would  otherwise  be  employed;  because  it  is  microcomputer-based, 
it  should  be  simpler  to  use,  and  because  the  mooring  legs  are  modeled 
using  catenary  formulas,  convergence  to  the  equilibrium  solution  should 
be  more  reliable. 

The  solution  algorithm  coded  in  SOUPLE  is  the  Newton-Raphson 
method.  Significantly,  the  algorithm  has  been  modified  to  enhance 
convergence  performance.  The  modification  consists  of  an  addition  of 
artificial  stiffness  to  the  mooring  system's  structural  stiffness 
matrix.  The  amount  of  artificial  stiffness  added  is  variable  in  two 
respects.  First,  the  user  is  required  to  make  decisions  on  the  appro¬ 
priate  amount  of  artificial  stiffness  during  the  input  data  phase  based 
on  his/her  experience  in  the  use  of  the  program  to  solve  more  difficult 
and  highly  nonlinear  problems.  For  more  well-behaved  problems  a  default 
value  may  be  selected.  Second,  the  amount  of  artificial  stiffness  dimin¬ 
ishes  automatically  after  each  iteration  of  the  solution  procedure  so 
that  only  a  minimum  amount  of  artificial  stiffness  exists  when  conver¬ 
gence  is  achieved. 

SOUPLE  does  have  its  limitations  in  static  problems,  and  it  looses 
some  of  its  formulation  elegance  when  handling  cable  inertia  forces  in 
dynamic  problems.  In  static  problems,  the  catenary  elements  used  to 
model  a  mooring  leg  must  each  remain  in  their  own  plane,  so  if  the  ver¬ 
tical  current  profile  is  not  coplanar,  the  mooring  leg  must  be  modeled 
as  a  series  of  smaller  elements  each  deforming  as  a  catenary  in  its  own 
plane.  Even  when  the  vertical  profile  is  coplanar,  but  not  uniform,  the 
user  must  anticipate  the  use  of  a  series  of  elements  to  model  the  varia¬ 
tion  of  drag  forces  with  depth.  This  can  be  done  easily  enough  by  spec¬ 
ifying  the  load  variation  on  the  mooring  leg  in  its  inivial  configura¬ 
tion.  However,  if  the  mooring  leg's  vertical  position  changes  during 
displacement  from  the  initial  position  to  the  final  position,  SOUPLE 
will  not  accomodate  the  changes  because  the  load  on  each  element  will 
remain  constant  irrespective  of  the  changing  depth  of  the  element. 

The  situation  in  dynamic  problems  is  similar  because  the  inertia 
forces  on  a  mooring  leg  will  vary  spatially.  Thus,  when  using  SOUPLE  to 
solve  dynamic  problems,  a  significantly  larger  number  of  catenary  elements 
will  be  needed  to  model  the  variation  in  inertia  forces  with  fidelity. 

This  mitigates  the  advantage  that  SOUPLE  would  normally  have  over  SEADYN* 
and  SEASTAR  due  to  its  super  catenary  element  formulation  in  contrast  to 
a  traditional  finite  element  formulation  for  modeling  compliant  struc¬ 
tures  . 


*While  SEADYN  also  has  this  catenary  element,  its  architecture  s  not 
designed  around  it.  SOUPLE  has  a  smaller  architecture  which  permits 
a  microcomputer  implementation. 
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